Dipolar Bosons in Triangular Optical Lattices: 
Quantum Phase Transitions and Anomalous Hysteresis 
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We study phase transitions and hysteresis in a system of dipolar bosons loaded into triangular 
optical lattices at zero temperature. We find that the quantum melting transition from supersolid 
to superfluid phase is first-order, in contrast with the previous report. We also find that due to 
strong quantum fluctuations the supersolid (or solid)-superfluid transition can exhibit an anoma- 
lous hysteretic behavior, in which the curve of density versus chemical potential does not form a 
standard loop structure. Furthermore, we show that the transition occurs unidirectionally along the 
anomalous hysteresis curve. 
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Ultracold atomic and molecular gases provide very 
clean and tunable systems to study various phenom- 
ena in condensed matter physics. Due to the remark- 
able control of physical parameters, one can simulate the 
physics of quantum many-body systems in regimes in- 
accessible to solid-state materials. Recently, two impor- 
tant developments have taken place in this area. First, 
experimental techniques for the preparation of ultracold 
gases with strong dipole-dipole interactions have been 
rapidly advancing over the last few years. This has been 
demonstrated by the realization of Bose-Einstein conden- 
sation (BEC) of ^^Cr atoms which have large magnetic 
dipole moments [l|, Q and by the creation of heteronu- 
clear (dipolar) polar molecules Secondly, trian- 

gular optical lattices of '^^Rb have been created exper- 
imentally, where the superfluid (SF) to Mott insulator 
transition has been observed @ . 

Stimulated by these experimental developments, we fo- 
cus on a system of ultracold dipolar bosons loaded into a 
triangular optical lattice. Due to its long-range nature, 
the dipole-dipole interaction coupled with the geometry 
of the triangular lattice can produce strong frustration. 
This setup provides an ideal venue for studying the in- 
terplay between strong frustration and quantum fluctu- 
ations. The studies of frustration have been carried out 
mainly in the field of magnetic materials. The frustration 
of spins can lead to exotic low-tcmperature spin states, 
such as spin glass, spin liquid, and spin ice (7H9|. 

It is well known that a system of lattice bosons with 
finite-ranged repulsion can be mapped, in the hard-core 
limit, onto a quantum spin-1/2 system with an XXZ- 
type anisotropy and a longitudinal magnetic field jioj . 
However, unlike usual magnetic materials, the exchange 
interactions of the mapped spin Hamiltonian are ferro- 
magnetic for the X and y components, which means that 



the frustration arises only from the coupling between the 
z components. Therefore the studies on strongly inter- 
acting bosons on frustrated lattices have great potential 
to pioneer new and intriguing phenomena not found in 
the regime of real spin systems and to provide a deeper 
understanding of geometrical frustration from a new per- 
spective. Thus, in this paper, we report the effects of 
quantum fluctuations and geometrical frustration in tri- 
angular optical lattices of dipolar bosons leading to quan- 
tum melting of supersolid (SS) or solid phases into SFs 
and to an anomalous unidirectional hysteresis in the den- 
sity versus chemical potential phase diagram. 

To capture the physics described above, we model 
dipolar bosons, for a sufficiently strong on-site interac- 
tion, by the following hard-core dipolar Bose-Hubbard 
model on the triangular lattice 
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J is the hopping amplitude be- 
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twecn nearest-neighbor sites, and /i is the chemical po- 
tential. We assume that the dipole moments are polar- 
ized by the external field in the direction perpendicular 
to the lattice plane. In this case, the interaction between 
the dipoles is isotroj)ic and can be well approximated by 
Vji = VS I \vj — r;| . Here, d is the lattice spacing. 

To study the quantum melting transitions and the ac- 
companying hysteresis curves, we use a cluster mean- 
field (CMF) method p^ - fisj . in which we can easily get 
the stationary points of the free energy not only for the 
globally stable state but also for metastable and unsta- 
ble states. To get meaningful results, one has to em- 
ploy a sufficiently large cluster for estimating the val- 
ues of mean fields. Although the use of the three-site 
cluster is convenient for tripartite lattices, the size 
iVc = 3 is too small to see clearly the effects of str ong 
quantum fluctuations. The work described in Ref. [1^ 
failed to capture some important features of the phase 
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FIG. 1: (color online) (a) The CMF-10 result for the 
ground-state phase diagram of hard-core bosons with nearest- 
neighbor repulsion on a triangular lattice in the [J /V, iJ,/V)- 
plane, where second and first-order phase transitions are in- 
dicated by thin blue and thick red lines, respectively, (b) 
Schematic pictures of the symmetries of p = 2/3 solid and 
upper SS states (left) and of p = 1/3 solid and lower SS 
states (right). The lattice sites are separated into triangular 
(A) and honeycomb (B) sublattices. (c) The filling factor p 
as a function of p/V for J /V = 0.1. (d) Cluster-size scaling 
(A 1) of the CMF data for the jump Ai in p at pi/V = 3. 
The linear fits of the three points [Nc ~ 6, 10, and 15) are 
fairly good. 



diagram, e.g., the existence of the direct solid-SF tran- 
sition at 1/3 < p < 2/3 and the discontinuity of the 
SS-SF transition which will be mentioned later. There- 
fore, using much larger-size clusters of triangular shape 
(up to Nc = 15), here we perform more complex, but 
more reliable, calculations |15|. 

To develop some intuition, we consider first hard-core 
bosons with only nearest neighbor interactions, i.e., we 
set Vji = V for nearest-neighbor bonds, and Vji = oth- 
erwise. The ground-state phase diagram of this simplified 
model contains a wide region of SS phase, in which long- 
range solid (crystalline) order and superfiuidity coexist, 
as well as the standard SF and solid phases [l6l - fl8| . Fig- 
ure Ufa) shows the ground-state phase diagram obtained 
by the ten-site CMF calculation (CMF-10). The SF state 
is characterized by the order parameter = '^j{aj)/M, 
where M denotes the number of lattice sites. The 
solid states with filling factors p = 1/3 and p = 2/3 
have the two-subblatice structures depicted in Fig.[TJb), 
which are characterized by pq = '^j{nj) cxp(iQ • Vj)/M 
with Q = (47r/3(i, 0). The filling factor is given by 
p = "^jifij) /M . In the SS states, |^| and |pq| have non- 
zero values simultaneously. We determined the boundary 
lines of first-order transitions from the Maxwell construc- 
tion in (J/y, x)-plane, where x = S(j,i>(°]"' +«3^«j)/^^- 

In Fig. [Tfc) we show the CMF-10 result for the fiUing 
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FIG. 2: (color online) (a) The J/V dependence of |pq|^ at 
p/V = 3.4 and p/V = 3.9. (b) The jump A2 in |pqp along 
the SS-SF (or solid-SF) transition line as a function of p/V . 



factor p as a function of p/V along the line of J/V — 0.1. 
At the particle-hole symmetry point p/V = 3, we can ob- 
viously see a finite jump, Ai ss 0.1 for J/V — 0.1, in the 
density curve, whereas previous numerical studies have 
shown that the density deviation |p— 1/2| at p/V = 3 is 
extremely small or even undetectable [13, [13 ■ To exam- 
ine this quantitative difference, we perform the infinite- 
size extrapolation of the results for different-size clusters 
[sec Fig. [IJd)] with the scaling parameter A defined by 
-^B / (^c X 2;/2) , where N-q is the number of bonds within 
the cluster and z = 6 (see Ref. [lH). The jump Ai de- 
creases with cluster size as expected. For example for 
J/V — 0.1, its value is strongly reduced from Ai w 0.1 
of CMF-10 to Ai « 0.01 in the limit A 1. The ex- 
trapolated value vanishes at J/V = 0.108, which means 
that the location of the triple point of the two SS and SF 
phases is shifted to {J/V, p/V) = (0.108,3) in 

The features of the phase diagram in Fig. [Ija) are in 
excellent agreement with those of QMC calculations [l^ ■ 
However, there is one major qualitative difference; our 
CMF results show that the transition between SF and 
SS is first-order (discontinuous). To confirm this, we 
plot IpqP as a function of J/V , calculated by CMF-10 
along the horizontal lines p/V = 3.4 and p/V = 3.9, 
in Fig. HJa). For the both values of ^/V , there is a fi- 
nite discontinuous jump A2 at the SS-SF transition, al- 
though that of the former is very small. The magnitude 
of jump decreases monotonically as the value oi p/V ap- 
proaches the particle-hole symmetry point. It vanishes 
at p/V = 3. As shown in Fig.[5Jb), the curve shows very 
little change with increasing cluster size, and thus we 
conclude that the SS-SF transition is first-order except 
for the critical point pi/V = 3. In contrast, this tran- 
sition seems to be continuous in QMC simulations, see 
Fig. 8 of Ref. [l^, where the authors concluded that the 
SS-SF transition is second- order. This discrepancy may 
come from a finite-size effect of QMC calculations, since 
the magnitude of the jump A2 near p/V = 3 is quite 
small as shown in Fig. [D^a). For values of p/V farther 
away from the particle-hole symmetry line p/V = 3, the 
discontinuous behavior is more pronounced, as shown in 
Fig-IHb), and can now be observed also within QMC [l9| . 

Let us discuss the hysteresis in the cycle of decreas- 
ing and increasing the chemical potential p/V . The sys- 
tem exhibits different hysteretic behaviors in three differ- 
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FIG. 3: (color online) (a) Magnified view of the low-density 
regions of Fig. [2a). We plot the limits of metastability of the 
SF phase (dashed green line) and the SS or solid phase (dash- 
dotted green line) in addition to the phase transition lines 
(thin blue and thick red lines) already shown in Fig.[TJa). (b) 
The usual hysteresis-loop and (c) anomalous unidirectional 
hysteresis behaviors. The thick solid, thin solid, and dashed 
lines represent ground, metastable, and unstable states, re- 
spectively. The inset in panel (c) shows the average density 
on each sublattice denoted by pA and pe, which satisfy the 
relation p — (pA + 2pb)/3. The corresponding points in the 
figures are marked with the same letters. 

ent ranges of J/V, which are defined by the thresholds 
{J/V)ci « 0.182, (J/y)c2 « 0.186, and (J/y)c3 « 0.220 
[marked by the dashed vertical lines in Fig. El^a)]. In 
the first region, J/V < {J/V)ci, accompanying the SF- 
solid transition, a typical hysteresis loop is formed in the 
{fi/V, p)-p\a,ne as indicated by the arrows in Fig [HJb). 
This is simply analogous to a conventional liquid-solid 
transition. In the second region, {J/V)ci < J/V < 
{J/V)c2, another hysteresis loop is formed around the 
SS-SF first-order transition point in addition to the loop 
around the SF-solid transition. 

Of particular interest is the third region, {J/V)c2 < 
J/V < {J/V)c3, in which the hysteresis exhibits 



an anomalous behavior. As an example, we show 
in Fig. ^c) the solution curves of the CMF-10 self- 
consistent equation in the {p/V, p)-p\a,ne for J/V = 
0.205. There are two first-order transitions, namely, be- 
tween the solid (at point e) and SF (f) states and be- 
tween the SS (h) and SF (i) states. Although the solu- 
tion branches corresponding to metastable SF and un- 
stable SS states apparently cross, the two states at the 
intersection arc not identical. This is clearly seen in the 
inset of Fig.[3]^c), where we plot the sublattice densities. 
Thus the solution curves are completely separated into 
the line of SF solutions and the twisted closed curve con- 
sisting of solid and SS solutions. In this case, we have 
an irreversible quantum melting transition, and once the 
solid order is melted at T = 0, it will remain melted. In 
this regime, the solid phases can be reached again only 
through thermal cycling. 

To illustrate this, let us assume that the system is ini- 
tially in a stable solid state located between points e and 
g in Fig.[3]Jc). When decreasing fi/V, although a SF state 
becomes energetically favorable below point e, the solid 
state remains metastable until it reaches point x. Below 
point X, the system is destabilized into the true ground 
state (namely the SF state). If wc increase n/V from the 
solid state, the system first undergoes a continuous tran- 
sition to the SS phase at point g, and then the metastable 
SS state is also destabilized into the SF state at point y. 
On the other hand, the situation drastically changes if 
we start from an initial state in the SF phase. We see 
in Fig. [3]Jc) that the globally stable SF solutions at low 
and high fi/V arc connected by the line of metastable 
SF solutions, which means that the SF state is stable for 
any p/V. Therefore, when wc decrease or increase fi/V 
starting from a SF state, the system remains in the SF 
phase even if the value of p,/V enters the region where a 
SS or solid state has the lowest energy. Thus the transi- 
tion in varying fj,/V occurs only from the SS (or solid) to 
SF phase, and the hysteresis curve does not form a stan- 
dard loop structure. Crucially, this behavior gets more 
pronounced and the region {J/V)c2 < J/V < {J/V)c.3 
becomes wider as the cluster size Nq increases, which 
means the robustness of the anomalous hysteresis when 
1. 

We also discuss the influence of the full long-range 
dipole-dipole interactions [l^l on our results. First, we 
perform a simple mean-field (MF) analysis for Eq. ([1]) 
with Vji ~ V(P/ \rj — r/| . As in the case of the nearest- 
neighbor interaction model, the MF phase diagram de- 
picted in Fig. IHa) includes large regions of p = 1/3 
and p = 2/3 solids and the two-sublatticc SS phase 
(named SSI) located between them. This is consistent 
with recent QMC calculations [ll|. Furthermore, the 
two-sublattice SS phase is stabilized also in the regions 
p < 1/3 and p > 2/3 (SSI'), and we find additional SS 
phases (SS2 and SS3) and solid phases with p = 1/2, 
p = 1/4, and p = 3/4 within the parameter range of 
Fig-SJa). To represent these phases, wc allowed for three- 
and four-sublatticc structures in minimizing the MF en- 



4 




JIV JIV 



FIG. 4: (color online) Ground-state phase diagram of hard- 
core bosons with full long-range dipole-dipole interactions on 
a triangular lattice in the ( J/y, /i/V)-plane, obtained from 
(a) MF theory and (b) CMF-10 method. Second and first- 
order phase transitions are indicated by thin blue and thick 
red lines, respectively. The dashed green lines represent the 
limit of metastability of the SF phase. The SSI and SSI' 
phases have the same sublattice structure as the SS phase in 
Fig. [Hb). (I-IV) Sublattice structures of the p = 3/4 solid 
and the nearby SS3 (I), of the p = 1/2 solid (II), of the SS2 
(III), and of the p = 1/4 solid and the nearby SS3 (IV). 

ergy. Next, using the CMF-10 method we take into ac- 
count quantum fluctuations around the MF solutions [sec 
FigHlJb)]. As in the case of nearest-neighbor interactions 
[Fig[IJb)], a noteworthy consequence of quantum fluctua- 
tions is that the SSl-SF boundary forms a dip around the 
particle-hole symmetry line. Moreover, a similar anoma- 
lous hysteresis described in the nearest-neighbor case also 
emerges here, and it is associated with the presence of the 
dip when (J/y)c2 < J/V < {J/VU [{J/V)c2 ~ 0.130 
and {J/V)c3 ~ 0.156]. 

Actual experiments of ultracold gases are performed in 
the presence of a trap potential, e.g., Vt(r) = ma;^|rp/2. 
Within the local-density approximation (LDA) , the effec- 
tive local chemical potential is given by jlj = fi — Vt{rj). 
We suggest that the anomalous hysteretie behavior can 



be confirmed by controlling (decreasing and increasing) 
jlj at the trap center via manipulation of, e.g., the fre- 
quency uj of the harmonic trap confining the dipolar 
gases. However, further analyses beyond LDA are still 
required to fully understand the behavior of the anoma- 
lous hysteresis in a trapped system. 

In summary, we have studied phase transition phe- 
nomena in a system of dipolar Bose gases loaded into 
a triangular optical lattice. Using a CMF method, we 
have found that the first-order transition between the SS 
(or solid) and SF phases can exhibit an anomalous hys- 
teretie behavior: in varying the chemical potential, the 
standard hysteresis loop structure does not appear, and 
the phase transition occurs only from the SS (or solid) to 
SF state. This unidirectional character is not predicted 
within the MF (classical) approach, since the boun dary 
of the SS-SF transition is given by a straight line jla |. 
Moreover, previous studies on a similar hard-core bo- 
son model with nearest-neighbor interactions for a square 
lattice have given only a standard hysteresis-loop behav- 
ior [2l|. Thus, the anomalous feature of the hysteresis in 
this system is attributed to the interplay between quan- 
tum fluctuations and the competition of interactions due 
to the frustrated geometry of the triangular lattice. 

More specifically, the most important point for the 
anomalous hysteresis behavior to emerge is the existence 
of the first-order SF-solid-(SS-)SF transition under vary- 
ing fJ,/V, as can be seen from Figs, ^a) and (c). Thus 
we expect to find analogous anomalous hysteresis in a 
wide range of systems exhibiting re-entrant first-order 
transitions (unless the system has no special symmetry 
pointjje.g., the Heisenberg point in the square-lattice 
case [2l|). Recently, it has been suggested that a spin 
supersolid state can be realized in some magnetic sys- 
tems [13 . Thus, by analogy, we also expect that such 
quantum spin systems will exhibit a similar anomalous 
hysteresis behavior as a function of magnetic field. 

We thank Grant-in-Aid from JSPS (D.Y.), KAKENHI 
(22840051) from JSPS (I.D.) and ARO (W911NF-G9-1- 
0220) (C.S.d.M.) for support. 
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A. Supplementary Material for "Dipolar 
Bosons in Triangular Optical Lattices: 
Quantum Phase Transitions and 
Anomalous Hysteresis" 

B. A. A large-size CMF method and cluster-size 
scaling 

The CMF method is convenient and efFective in under- 
standing the physics of ordered states including metasta- 
bihty phenomena since all stationary points of the free 
energy can be obtained not only for the globally stable 
solution. Moreover, this method is free from "finite-size 
effects" and "error bars." However, in order to get rea- 
sonable results, one has to treat sufficiently large clusters 
as a reference system, especially when strong fluctuations 
exist in the system. In this work, we have used a series of 
clusters of Nq = 3, 6, 10, and 15 sites given in Table [Sll 

For a large cluster, we cannot calculate the value of 
free energy of the system directly from the CMF for- 
malism, unlike the standard MF theory. This is a im- 
portant problem for the system considered in this work, 
which exhibits first-order phase transitions in a wide re- 



TABLE SL A series of clusters used in our calculations. The 
values of A^c, A^b, and A are also listed. 

MF CMF-3 CMF-6 



A A 



A^c 1 3 6 

Nb 3 9 

A 1/3 1/2 



CMF-10 CMF-15 Exact 




Ac 10 15 oo 
Nb 18 30 oo 
A 3/5 2/3 1 



gion of parameters. To overcome this problem, we have 
used the Maxwell construction in {J/V, x)-plane to de- 
termine the phase boundary. The quantity x is defined 
by X = -I- a\aj)/M, and the energy differ- 

ence between the states at Jq and Ji (when /x and V are 
fixed) is given by — /j ^ x{J)dJ- Note that the Maxwell 
construction in (/i/V, /9)-plane [l| is inapplicable in this 
case because of, indeed, the existence of the anomalous 
hysteresis; as shown in Fig. 3(c), the solution curves are 
separated into two groups for {J/V)c2 < J/V < {J/V)c3, 
and thus one cannot estimate the energy difference be- 
tween states of different groups by using the integration 
of the density p over the chemical potential fJ-/V. 

Moreover, to estimate the expectation values (•••), one 
should take the average not only over all internal sites 
within the cluster, but also all possible choices of how 
to embed the cluster itself in the background sublat- 
tice structure. For example, when we assume the two- 
sublattice -s/S x -\/3 ordering in the ten-site CMF ap- 
proximation (CMF-10), we have three choices of clusters 
shown in Fig.lSll While we treat exactly the interactions 
within the cluster, the interactions between the cluster 
and the rest of the system are also included via effective 
fields acting at the cluster edge. The effective fields are 
determined self-consistently via the expectation values of 
the operators (aj) and {nj), in a standard manner 0. 

In the main text, we have performed a cluster-size scal- 
ing of the CMF results for the amplitude of a discontinu- 
ous jump Al at the transition between the two SS phases 
[see. Fig. 1(d)]. To extrapolate the value of Ai, we intro- 
duced the scaling parameter A defined by A'^b / (-^c x ^/2) , 
which varies from to 1. Here, A^b is the number of bonds 
within the cluster and z is the coordination number for 
the lattice, i.e., z = 6 in this case. The denominator 
means the number of bonds of the original lattice per A^c 
sites, and hence the parameter A provides an indication 
of how much the correlation effects between particles are 
taken into account in the cluster. The value of A for each 
cluster is shown in Table [Sll In Fig. 1(d) we performed a 
linear extrapolation toward A = 1 using the three points 
of CMF-6, -10, and -15, since the three-site cluster is too 
small for scaling. The linear fits of these three points are 
fairly good, as seen in Fig. 1(d). At J/V ^ 0.108, the 




FIG. SI: Schematic picture of the ten-site clusters embedded 
in the background two-sublattice x \/3 structure. 



6 



p=2/3 Solid 

Vi (b) »' 








SSA 






SSB 





p= 1/3 Solid i^'^ 

(c) "r 



^ ...I 













/ SSB 


ssa\ 







SSB 



SSC 



SSA 



FIG. S2: The /x/V dependences of (a) the average density 
p, (b) the sohd order parameter |pq|, and (c) the SF order 
parameter |*| for J /V = 0.1, obtained by the CMF-10 cal- 
culations, (d) Magnified view of the region near ^/V = 3 in 
(a). The thick solid, thin sohd, and dashed lines represent 
the solutions of the ground, metastable, and unstable states, 
respectively. The lower panels show the sublattice structures 
of SSA, SSB, and SSC. 



jump Ai vanishes in the limit A — > 1, which means that 
the triple point of the two SS phases and the SF phase 
is strongly reduced from the MF value J/V = 0.25 Q 
to J/V — 0.108 due to the quantum fluctuations. This 
estimate for the triple point is in fairly good agreement 
with the recent QMC prediction Q. 



called SSA for distinction) and the low-density SS phase 
(SSB), since the main focus was the hysteresis occurring 
near the transition between solid (or SS) and SF phases. 
We present here some additional results for the properties 
of the SSA-SSB transition. 



In Figs. EKa-d), the results of the CMF-10 for the 
filling factor p, the solid order parameter |pq| with Q = 
(47r/3(i, 0), and the SF order parameter |\['| are shown as 
functions of fi/V for J/V = 0.1. As shown in Fig.[S2ja), 
the ground-state density curves of SSA and SSB are not 
connected at the center, while the two order parameters 
in Figs. IS2rb') and (c) do not show such discontinuity. 
The latter is because of the particle-hole symmetry at 
fi/V = 3. The discontinuity in the density indicates the 
first-order character of the transition between the two SS 
phases, although the amplitude of the jump is strongly 
reduced in the limit A — >■ 1 [see, Fig. 1(d)]. 



In addition, we plot in Fig. IS2rd') the solution curve 
of the CMF-10 equation, including the metastable and 
unstable solutions, near the SSA-SSB phase boundary. 
The solution lines of SSA and SSB are extended due 
to the metastability, and connected by the line of the 
three-sublattice SS state, called SSC We found 

that unstable SSC solutions (corresponding to unstable 
stationary points of the free energy) exist in the vicin- 
ity of half filling, although there is no region of the SSC 
phase in the ground-state phase diagram. The density 
versus chemical potential curve has negative curvature 
(negative compressibility) in the SSC state, which indi- 
cates the emergence of phase separation at a fixed density 
near p = 1/2; the system separates into a mixture of SSA 
and SSB phases. 



C. B. Additional data for the transition between 
two SS phases 

In the main text, we made only a passing reference to 
the transition between the high-density SS phase (often 



Note that in the case of infinite-range dipole-dipole in- 
teractions, there is a finite region of stable SSC phase in 
the phase diagram (marked SS2 in Fig. 4), although it is 
remarkably narrow for CMF-10. It is an interesting fu- 
ture issue to determine whether the SSC phase is actually 
stable or not for the infinite-range interaction model. 
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